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ABSTRACT 

We model the cosmological co-evolution of galaxies and their central supermassive black 
holes (BHs) within a semi-analytical framewor k developed on the o utpu ts of the Millennium 
Simul ation. This model, described in detail in Croton et alJ (l2006l) and iDe Lucia & Blaizo3 
(12001) . introduces a 'radio mode' feedback from Active Galactic Nuclei (AGN) at the centre 
of X-ray emitting atmospheres in galaxy groups and clusters. Thanks to this mechanism, the 
model can simultaneously explain: (i) the low observed mass drop-out rate in cooling flows; 
(ii) the exponential cut-off in the bright end of the galaxy luminosity function; and (iii) the 
bulge-dominated morphologies and old stellar ages of the most massive galaxies in clusters. 
This paper is the first of a series in which we investigate how well this model can also re- 
produce the physical properties of BHs and AGN. Here we analyze the scaling relations, the 
fundamental plane and the mass function of BHs, and compare them with the most recent 
observational data. Moreover, we extend the semi-analytic model to follow the evolution of 
the BH mass accretion and its conversion into radiation, and compare the derived AGN bolo- 
metric luminosity function with the observed one. While we find for the most part a very good 
agreement between predicted and observed BH properties, the semi-analytic model underes- 
timates the number density of luminous AGN at high redshifts, independently of the adopted 
Eddington factor and accretion efficiency. However, an agreement with the observations is 
possible within the framework of our model, provided it is assumed that the cold gas fraction 
accreted by BHs at high redshifts is larger than at low redshifts. 



1 INTRODUCTION 



Over the last years, several observations have demonstrated that su- 
permassive black holes likely reside at the centres of all spheroidal 
galax ies (see e.g. iKormendv & Richstond 1 19951 ; IRichstone et al.l 
1 19981) . Even more interestingly, their prope rties appear to strongly 
correlate with those of the ir host gal axies (Magorrian et al. 199o; 
Ferrarese & Merritj |2000|; lOebhardt et aL .2 000 ; Gra ham et al 



20011 ; ITremaine et ail I2OO2I; iMcLure & Dunlop 2002; Baes et all 
2003"; 'Marconi et al.' '2004'; Harin g & Rixl 120041 ; IFeoli & Meld 
2007; Graham & Driver 2007) and, apparently, a lso with the ones 
of the whole host dark matter (DM) haloes ( lFerrares3 1 20021 ; 
iBaes etal.|[2003l ; iFerrarese & Fordl2005l) . Although it is not yet 
clear which of th ese relations is "more fundamental" (see e.g. 
iNovak et alj |200^. they reasonably suggest a close link between 
the assembly history of t he BHs and the cosmo logical evolution 
of galaxies. Most recently, iHopkins et id] ( 12007 j) have shown that 
these relationships are not independent and could be interpreted 
as different projections of a BH fundamental plane, analogous to 
the fundamental plane of elliptical galaxies. The striking similarity 



between these two fundamental planes is another clue that galaxy 
spheroids and BHs do not form and evolve independently. 

The paradigm that AGN are powered by mass accretion onto 
BHs (Salpeter 1964; Lvnden-Bell 1969) has also received very 
strong support from spectroscopic and photometric observations 
of the stellar and gas dynamics in the central regions of local 
spheroidal galaxies and bulges. Moreover, by estimating the total 
energy radiated by AGN during their whole life, it can be shown 
that nearly all the mass in B Hs has been accumulated during peri- 
ods of bright AGN activity jSoltanlll982l) . implying that the com- 
mon physical process which produces galaxy spheroids and BHs 
must also be responsible for triggering bright AGN. 

Such a cosmological co-evolution of BHs, AGN and 
galaxies is expected in the standard framework, in which 
cosmic structures grow hierarchically via gravitational insta- 
bility and merging events destabilize the gas at the galaxy 
centres, triggering star formation and BH mass accretion. In 
order to investigate this complex scenario, several models have 
been developed, based on either pure analytic approximations 
(see, e.g., Efstathiou & Rees 1988; Haehnelt & Rees 199f 
iHaiman&Loebll 19981 ; IPercival & Milleill 19991 ; IHaiman & Meno 
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to the availability of unprecedented 
computational powe r, fully numerical mo dels have also become 
available (see e.g.. loi Matteo et alj 12005: Springel et all l2005l ; 
iLi et al.ll2007l ; ISiiacki et alj2007l ; lDi Matteo et alj|2007h . 

Simple analytic models in which AGN activity is only trig- 
gered by DM halo major mergers succeeded in quantitatively de- 
scribing the observed evolution of the AGN number counts and lu- 
minosity at all but low redshifts, provided that some mechanism is 
advocated to inhibit accretion within massive haloes hosting bright 
AGN. However they fail in reproducing the observed AGN clus- 
tering at high redshifts (Marulli et al. 2006). Slightly more sophis- 
ticated semi-analytic models in which the halo merger history and 
associated BHs are followed by Monte Carlo realizations of the 
merger hierarchy, while the baryonic physics is neglected as well, 
can correctly reproduce both the AGN luminosity and clustering 
function at z > 1 (Marulli et al. 2006), but the number density of 
faint AGN is significantly below observations, a clear indication 
that DM halo mergers cannot constitu te the only trigger to accretion 
episodes in the local BH population dMarulli et alj|2007l) . and that 
in order to properly describe the cosmological evolution of BHs and 
AGN, the main baryonic phenomena involving the gas contents of 
DM halos cannot be neglected. 

This complication is reminiscent of the one found in the de- 
scription of galaxies, where the well-known mismatch in shape be- 
tween the predicted distribution of DM halo masses and the ob- 
served distribution of galaxy luminosities requires the considera- 
tion of complex baryonic phenomena like, for instance, cooling in- 
efficiencies to reduce gas condensation in massive structures, or 
supernova ( White & Rees 1978; White & Frenk 1991) and stellar 
kinetic feedback dFontanot et al.ll2006l) to remove cold gas in low 
mass systems, as we ll as photoionisat ion heating to suppress the 
formation of dwarfs ( lEfstathiou|[l992h . Cooling effects alone are 
however too weak to produce the bright end cut-off of the lumi- 
nosity function, and it seems to be mandato ry to include addi- 
tional feedback processes in massive halos (e.g.lBenson et al .I2OO3I ; 
lFontanotetai]|2006l ; ICroton et al.ll2006l ; Iciotti & OstrikeJl2007h . 
Standard models of galaxy formation face two additional prob- 
lems; i) the persistence of a hot gas atmosphere at the centre of 
most galaxy clusters despite the fact that the loc al cooling time is 
much shorter than the age of the system (see, e.g.lCowie & Binnevj 
19771; lFabian&NulseJ 'l977'; 'Peterson et al.' '2OOI'; iTamura et al.' 



2001 



iFabian etall2 003; McNamara et al. 2005; Morandi & Ettori 



2007 *. and references therein), and ii) the fact that most massive 



galaxies, typically ellipticals in clusters, are made of the oldest stars 
and so finis hed their star formation earlier than lower mass galax- 
ies (s ee, e.g. lCowie et alJll996l ; lDe Lucia et al.ll2006l ; [cimatti et all 
i2ooa and references therein). 

In this paper we study the cosmological co-evolution of 
galaxies and their central BHs using a semi-analytical model 
developed on the o utputs of th e Mill enni um Simulation and 
described in detail in lCroton et al.l l l2006h and lDe Lucia &. BlaizotI 
( I2OO7I) . In this scenario, radio mode feedback from AGN at 
the centre of galaxy groups and clusters is invoked to prevent 
significant gas cooling in large halos, thus limiting the mass 
of the central galaxies and preventing them from forming stars 
at late times when their mass and morphology can still change 



through mergers. Thanks to this mechanism, ICroton et al. I l l2006l) 
demonstrated that such a model can simultaneously explain the 
low observed mass drop-out rate in cooling flows, the exponential 
cut-off in the bright-end of the galaxy luminosity function, and 
the bulge-dominated morphologies and stellar ages of the most 
massive galaxies in clusters. 

Here we are interested in investigating how well this model 
can also reproduce the statistical properties of BHs and AGN. 
To do that, we extend the original model by adding new semi- 
analytical prescriptions to describe the BH mass accretion rate in 
the accretion episodes triggered by galaxy mergers, which fuel the 
quasar mode, and their conversion into radiation. We then analyze 
the scaling relations, the fundamental plane and the mass function 
of BHs, and compare them with the most recent observational 
data available. Finally, we compare the predicted AGN bolometric 
luminosity function with the observed one, and propose some 
modifications to the original semi-analytic assumptions to better fit 
the data. 

The paper is organized as follows. In Section|2l we briefly de- 
scribe the main aspects of our semi-analytic model and illustrate the 
new equation introduced to describe the BH mass accretion in the 
quasar mode in more detail. In Section |3] we compare the model 
predictions with the best observational data available for the BH 
and AGN populations. Finally, in Section|4]we summarize our con- 
clusions. 



2 THE MODEL 

Our semi-analytic model for the co-evolution of DM haloes, galax- 
ies and their central BHs consists of three ingredients, that we de- 
scribe separately in this section: a numerical simulation to obtain 
the merger history of the DM haloes, a set of analytic prescriptions 
to trace the evolution of galaxies within their host haloes and a set 
of recipes to follow the BH accretion history and the AGN phe- 
nomenon. 



2.1 Numerical simulation 

In this work we use the outputs of the Millennium Simulation, 
which followed the dynamical evolution of 2160-' ~ lO''* DM par- 
ticles with mass 8.6 x 10^ /j^'M0 in a periodic box of 500/;^ 'Mpc 
on a side, in a ACD M "concordance" cosmological framework 
dSpringel et alj200^ . The computational box is large enough to in- 
clude rare objects such as quasars or rich galaxy clusters, the largest 
of which contain about 3 million simulation particles at z = 0. At 
the same time, the mass resolution is high enough to resolve the 
DM halo of 0. IL* galaxies with ^ 100 particles. The short-range 
gravitational force law is softened on the co-moving scale 5/i~'kpc 
(Plummer-equivalent) which may be taken as the spatial resolution 
limit of the calculation. The cosmological parameters (the mat- 
ter density parameter fim = 0.25, the baryon density parameter 
fit, = 0.045, the Hubble parameter h = //o/100kms"'Mpc"' = 
0.73, the cosmological constant contribution to the density param- 
eter flyv = 0.75, the primordial spectral index « = 1, and the power 
spectrum normalization Gg = 0.9), are consistent with determina- 
tions from the combined analysis of the 2-de gree Field Galaxy Red- 
shift Survey (2dFGRS) dColless et al.ll2()0ll) and first-year W MAP 
data dSpergeletal]|2003l) , as shown bv lSanchez et alj d2006h . We 
recall that the more recent analysis of the WMAP 3-year data 
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dSpergel et alll2007h suggests slightly different values (in particu- 
lar smaller values for Sl^a, CSg and n). However, as demonstrated by 
IWang et alj t-QQj). due to the current modelling uncertainties, it is 
not possible to distinguish the two WMAP cosmologies on the ba- 
sis of the observed galaxy properties, since the variations induced 
by acceptable modifications of the free parameters of the galaxy 
formation model are at least as large as those produced by the vari- 
ation in the cosmological parameters. 

The Millennium Simul ation was carrie d out with a special ver- 
sion of the GADGET-2 code ( ISpringell2005h . optimized for very low 
memory consumption, at the Computing Centre of the Max-Planck 
Society in Garching, Germany. We make use of hierarchical merg- 
ing trees extracted from this simulation which encode the full for- 
mation history of DM haloes and subhalos, previously identified 
with, respectively, a friends-of-friends (FOF ) group-finder and a n 
extended version of the SUBFIND algorithm jSpringel et alj|200lh . 
These trees constitute the backbone of our semi-analytic model, 
which is implemented during the post-processing phase: this allows 
us to simulate the wide range of baryonic processes occurring dur- 
ing the formation and evolution of galaxies and their central BHs. 



2.2 Galaxy evolution 

We use t he galaxy fo rmation model of ICroton et al. I ( I2OO6I) as up- 
dated by IDe Lucia & Blaizot i.2QQl])- Although not in agreement 
wit h some properties of the red and blue galaxy popu lations (see, 
e.g- IWeinmann et al.l2006l : lKitzbichler & Whitell2007h . this model 
is able to reproduce the overall observed properties of galaxies, i.e. 
the relations between stellar mass, gas mass and metallicity, the lu- 
minosity, c olour and morphology distributions (Croton et al. 2006; 
iDe Lucia et al. 2006), the two-point galaxy correlation functions 
I Springel et al..i2005t) , and the global galaxy luminos ity and mass 
functions at high redshift jKitzbichler & Whiteir2007l) . We refer to 
the original papers for a full description of the numerical imple- 
mentation of the model. In the following, we briefly recall the treat- 
ment of the physical processes involved in the galaxy evolution, and 
describe the prescriptions for the BH growth and the AGN evolu- 
tion. 

Following the standard paradigm set out by IWhite & FrenkI 

(Il99lh and a d apted to high-resolution N-body simulations by 
ISpringel et all l l200lh . we assume that as a DM halo collapses, 
a fraction fj, = 0.17 of its mass is in the form of baryons and 
collapses with it, c onsistent with the first-year WMAP result 
dSpergel et ai]|2003l) . Initially, these baryons are in the form of a 
diffuse gas with primordial composition, but later they include gas 
in several phases as well as stars and heavy elements. Convention- 
ally, with the simplifying assumption of an ideal gas which cools 
isobarically, the cooling time of the gas is computed as the ratio of 
its specific thermal energy to the cooling rate per unit volume. 



distribution. 



^cool 



3 liinpkT 
2p^(r)A(r,Z) 



(1) 



where fAnip is the mean particle mass, k is the Boltzmann constant, 
Pe(r) is the hot gas density, and A(T,Z) is the cooling function 
dSutherland & Dopitall993l : lMaio et alj2007h . Equation {D is valid 
at temperature higher than ~ lO'* K, where hydrogen and helium 
remain ionized and the number of particles remains approximately 
constant. 

We assume the post-shock temperature of the infalling gas to 
be the virial temperature of the halo, T = 35.9 (Vyir/kms^' )^ K, 
where Vyir is the halo virial velocity. Moreover, we assume that 
the hot gas within a static atmosphere has a simple 'isothermal' 
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where mj,ot is the total hot gas mass associated with the halo and is 
assumed to extend to its virial radius ^^ii- 

In order to estimate an instantaneous cooling rate onto the 
central object of a halo, given its current hot gas content, we de- 
fine the cooling radius, r^Qoh as the radius at which the local cool- 
ing time (assuming the structure of equatio n (|2)) is equal to the 
halo dynamical time , j^yir/Vyii- = l lSpringel et"al]|200ll: 

iDe Lucia et alj|2004l : [Croton et aLll2006h : here H{z) represents the 
redshift evolution of the Hubble constant. The cooling rate can then 
be determined through the following continuity equation, 

"'cool = 43tp^(rcool)rcool''cool • (3) 

More details abou t our cooling prescriptions can be found in 
ICrotonetal] ( l2006l) . 

The photo-ionization heating of the intergalactic medium 
suppresses the c oncentration of baryons in shallow potentials 
(Efstathioui il992l) . and can be responsible of t he inefficient ac cre- 
tion and cooling in low-mass haloes. Following lGnedinl ( l200iDh . we 
model the effect of such photo-ionization heating by defining a 
characteristic mass scale, Mp, below which the gas fraction ft, is 
reduced relatively to the universal value: 



(4) 



[l+0.26MF(z)/Mvir]3 

We adopt the Mp{z) parameterization of iKravtsov et al.l ( |2004|) . 
which results in a filtering mass Mf of 4 X 10^ Mq at the reioniza- 
tion e poch, and 3 x lO'^My, by the present day (but see lHoeft et al.l 
l2006l) . 

In the semi-analytic framework we use in this work, the star 
formation is assumed to occur at a rate given by: 

= asF(mcold -'''crit)Adyn.disc > (5) 

where nicoid is the cold gas mass, fjyn.disc is the dynamical time 
of the galaxy, defined as the ratio between the disk radius and the 
virial velocit y, merit correspond s to a critical value for the gas sur - 
face density ( lKauffmamilll996l : lKennicutill998l : IMo et alj|l998h . 
and asF = 03 controls the efficiency of the transformation of cold 
gas into stars. Massive stars explode as supemovae shortly after star 
formation events and are assumed to reheat a gas mass proportional 
to the mass of stars: 

AnJreheated = EdiskAm*, (6) 

where we set the free parameter Cdisk = 3.5 based on the observa- 
tional data. The energy released by an event which forms a mass 
Am* in stars is assumed to be: 

A£sN = 0.5ehaioAra*v|v, (7) 

where 0.5^5]^) is the mean supernova energy injected per unit mass 
of newly formed stars, and Ehaio represents the efficiency with 
which this energy is able to convert cold interstellar medium into 
hot, diffuse halo gas. The amount of gas that leaves the DM halo 
in a "super-wind" is determined by computing whether excess SN 
energy is available to drive the flow after reheating of material to 
the halo virial temperature. 

We model the disk i nstabilities using the analytic stability cri- 
terion of IMo et al] ( Il998l) : the stellar disk of a galaxy becomes un- 
stable when the following inequality is met: 

7? ? M/2 ^ 1 • 

(Gmdisk/'-disk) ' 
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At each time-step we evaluate the left-hand side of equation ^ 
for each galaxy, and if it is smaller than unity we transfer enough 
stellar mass from disk to bulge (at fixed vq) to restore stability. 

In the Millennium Run, substructures are followed down to 
masses of 1.7 X IO'^/i^'Mq, so that we can properly follow the 
motion of galaxies inside their hosting DM haloes until tidal 
truncation and stripping disrupt their subhalos at this resolution 
limit. At this point, we estimate a survival time for the galax- 
ies u^ing_Jheir_current orbit and the dynamical friction formula 
of Birmev & Tremairi3 1 19871) multiplied by a factor of 2, as in 
lOe Lucia & Blaizot 1 20071) . After this time, the galaxy is assumed 
to merge onto the central galaxy of its own halo. Galaxy mergers 
induce starburst which we describe using the "colli sional starburst" 
prescription introduced bv lSomerville et al] ( 1200 ih . In this model, 
a fraction eturst of the combined cold gas from the two merging 
galaxies is turned into stars as follows: 

fiburst = Pburst(msat/'"central )"'"'"" , (9) 

where the two parameters are taken as (X\^am = 0-7 and pbuist = 
0.56, appro priate for merger mass ratios ranging from 1:10 to 1:1 
( ICo)3l2()04h . 

2.3 BH mass accretion and AGN 

2.3.1 The 'radio mode' 

When a static hot halo has formed around a galaxy, we assume that 
a fraction of the hot gas continuously accretes onto the central BH, 
cau sing a low-energy 'r adio' activity in the galaxy centre. Follow- 
ing ICroton et al.l J2OO 0), the BH mass accretion rate during these 
phases is postulated to scale as follows: 

/ MbH \ (hoi\ ( Vvir 

\IQ^Mq) VO.l ) V200kms-1 



^BH.R = KaGN 



(10) 



where Mbh is the BH mass, /j,ot is the fraction of the total halo mass 
in the form of hot gas, and k^gn is a free parameter set equal to 
7.5 X lO^^Moyr^' in order to reproduce the turnover at the bright 
end of the galaxy luminosity function. Since /hot is approximately 
constant for VVh It 150kms^', the dependence of otbh.r on this 
quantity has a little effect. Note that the accretion rate given by 
equation l IlOl l is typically orders-of-magnitude below the Eddington 
limit. In fact, the total mass growth of BHs in the radio relative to 
the quasar mode discussed below is negligible. 

It is also assumed that the radio mode feedback injects en- 
ergy efficiently into the surrounding medium, which can reduce or 
even stop the cooling flow in the halo centres. The mechanical heat- 
ing generated by this kind of BH mass accretion and described as 
^BH = EJ^bhc^, where e = 0.1 is the accretion efficiency and c is 
the speed of light, induces a modified infall rate of the following 
kind: 



"cool — ™cool - 
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For consistency we never allow rir'^^^i to fall below zero. In this sce- 
nario, the effectiveness of radio AGN in suppressing cooling flows 
is greatest at late times and for large values of the BH mass, which 
is required to successfully reproduce the luminosities, colours and 
clustering of low-redshift bright galaxies. 

2.3.2 The 'quasar mode' 

In our model BHs accrete mass after a galaxy merger both through 
coalescence with another BH and by accreting cold gas, the lat- 



ter being the dominant accretion mechanism. For simplicity, the 
BH coalescence is modelled as a direct sum of the progeni- 
tor masses, thus ignoring gravitational wave losses. Following 
iKauffmann & Hae hnelt (2000), we assume that the gas mass ac- 
creted during a merger is proportional to the total cold gas mass 
present, but with an efficiency which is lower for smaller mass sys- 
tems and for unequal mergers: 



AMbh.q ^ 



/bh '"cold 



l + (280kms-VVvir)2 



where 



/bh = /bh ('"sat/'Kcentral) 



(12) 



(13) 



and /bh ~ 0.03 is chosen to reproduce the observed local Mbh — 
Mbuige relation. Thus, any merger-induced perturbation to the gas 
disk (which might come from a bar instability or a merger-induced 
starburst) can in principle drive gas onto the central BH. However, 
the fractional contribution of minor mergers is typically quite small, 
so that accretion driven by major mergers is the dominant mode of 
BH growth in our scenario. This kind of accretion, which we call 
quasar mode, is also closely associated with starbursts, which occur 
concurrently. We do not model feedback from the quasar activity in 
the current model, but it can be approximately represented by an 
enhanced effective feedback efficiency for the supernovae associ- 
ated with the intense starburst. 



2.3.3 AGN luminosity 

The output of the m odel summarized h it herto, called DeLu- 
cia2006a catalogue dOe Lucia & Blaizo3 l2007h . is publicly 
available at http://www.mpa-garching.mpg.de/millenniumj 
dLemson & Virgo Consortium .2006,) . In this default model, 
for simplicity, the BH mass accretion triggered by each merger is 
implemented as an instantaneous event and the BH seed masses 
are set equal to zero. 

In order to study the evolution of AGN inside this 
cosmological framework, we improve the original model of 
iDe Lucia & Blaizo J toot by adding new semi-analytical pre- 
scriptions to describe the BH mass accretion rate during each 
merger event in the quasar mode, and its conversion into radiation. 
In this implementation, BHs do not accrete mass instantaneously. 
Instead, the accretion is coupled to the light curve model adopted. 
If a galaxy undergoes a merger while the central BH is still accret- 
ing mass from a previous merger, the cold gas still to be accreted is 
added to the new gas reservoir, and the accretion re-starts under the 
new physical conditions. In Sect. l3.1.4l we show that the BH scaling 
relations are weakly affected by this change. We use the following 
definitions to parameterize the bolometric luminosity emitted by 
accretion onto BHs, as a function of the accretion efficiency, e, and 
the Eddington factor, /Edd(0 := i'bol(0/iEdd(f)' 

ibol(f) = y^^Bh(Oc^ 

= ./Edd(O^Edd(0 =./Edd(f)— r^^f^, 

^Edd 



•rflnMBH(f) = 



dt 



(14) 



where L^^d is the Eddington luminosity, t^^i^ = ajc/{4TZmpG) 
0.45 Gyr and fef(f) = yj^jTy is the e-folding time (?ef = fsalpet. 
if/Edd= 1)- 



l-E /Edd(') 

No strong observational constraints are available for £ and 
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Figure 1. The time evolution of /edd (top panel), Mbh (central panel) and 
i-bol (bottom panel) for our three lightcurve models (/ (blue solid lines), // 
(red short-dashed lines) and /// (green lines)), for an illustrative case of a BH 
of mass Mbh = IO'M.j accreting a mass AMbh.q = 5 x 10* Mo, starting at 
z = 3. The three green curves, showing our model ///, have been obtained 
by setting f = 0.5 (short dashed), 0.7 (dotted-long dashed) and 0.9 (short 
dashed-long dashed). 



/Edd> the parameters that regulate the BHs powering the AGN 
and, more importantly, if and how they depend on redshift, BH 
masses, AGN luminosities and so on. However, some observa- 
tions a^_zj;-^_OJndicate that 0.04 < e < 0.16 and 0.1 < /edd < 
1.7 jMarconi et alj|2004l) . Furt hermore, it has been suggested that 
/prtri may depends on redshift jShankar et alj|20()4h and BH mass 
dNetz er & Trakhtenbrot 2007). In this paper, for si mplicity, we do 
not fo llow the evolution of the BH spins (see, e.g. IVolonteri et al.l 
I2OO7I and references therein) and we take a constant mean value 
for the accretion efficiency of e = (e) = 0.1 at all redshifts. 

For /Eddi which determines the lightcurves associated with in- 



dividual quasar events, we consider instead three different prescrip- 
tions: 

• I- /Edd = 1' the simplest possible assumption. Here the quasar 
is either 'on' at its maximum Eddington luminosity, or 'off'. 

• //: 



/Edd(z) 



z<3 



(15) 



/Edd.O 

/Edd,0-[(l+z)/4]l-4 

with /Edd.O = 0-3, as suggested bv lShankar et al.l ( l2004l) to match 
the BH mass function derived from a deconvolution of the AGN 
luminosity function and the local BH mass function. 

• ///: based on the analysis of self-consist e nt hyd rodynamical 
simulations of galaxy mergers. iHopkins et ai] ( I2OO5I) noticed that 
the light curves of active BHs are complex, showing periods of 
rapid accretion after "first passage" of the merging galaxies, fol- 
lowed by a long-lasting quiescent phase, then a transition to a 
highly luminous, peaked quasar phase, finally a fading away when 
quasar feedback expels gas from the remnant's centre in a self- 
regulated mechanism after the BH reaches a critical mass. In spite 
of this complexity, as a first order approximation, the typical evo- 
lution of an active BH can be simply described as a two-stage pro- 
cess of a rapid, Eddington-limited growth up to a peak BH mass, 
proceeded and followed by a much longer quiescent phase with 
lower Eddington ratios. In this latter phase, the average time spent 
by AGN per logarithmi c luminosity interval can be approximated 
as jHopkins et"al]|2005l) 



df 



(16) 



dlnLboi V10%x 

where fg = tqijJ > IO^Lq) and tQ{L' > L) is the total AGN 
lifetime above a given luminosity L; tg ~ 10^ yr over the range 
lO^Lfy, < Lh„i <is|ak- In the range lO^'La < Lpeak < lO'^Lg, 
iHopkins et al ] ( I2OO5I) found that a is a function of only the AGN 
luminosity at the peak of its activity, Lpeak^ given by a = —0.95 + 
O.321og(Lpeak/lO'^^0)> with a = —0.2 (the approximate slope of 
the Eddington-limited case) as an upper limit. We here interpret 
the Hopkins model as describing primarily the decline phase of the 
quasar activity, after the black hole has grown at the Eddington rate 
to a peak mass Mbjj peak = A^BH(fin) + J' - AMbh.q- (1-e), where 
^BH('in) is the initial BH mass and AMbh.q is the fraction of cold 
gas mass accreted. Here !f is an additional free parameter, in the 
range ^ T ^1. For T = 1 the BH emits at the Eddington rate. In 
the opposite limit (f = 0) the AGN reaches instantaneously a peak 
luminosity, and the whole light curve is described by equation 1 II6I 1. 
We found that f = 0.7 is the value that best matches the AGN lu- 
minosity function. We note that this interpretation of the Hopkins 
model is plausible but not unique, as part of the time described by 
equation l ll6t could also be associated with the rising part of the 
lightcurve. 

From equation 1 II6I 1 and with the following definition 



/Edd(0 



we can derive: 



d/Edd(0 
df 



-/Edd(0 



-'peak 



rl-a 
J Edd 



-'peak 



^peak 



fa 



I ^peak 



1/a 



(17) 



(18) 



(19) 



/Edd.O ' I lO'^L 

Here we neglected the absolute value of a present in equation l ll6t . 
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Figure 2. Starting from the upper left panel down to the bottom right one, scaling relations between the masses of the central BHs in the simulated galaxies 
with six different properties of their hosts: the K- and B-band bulge magnitude (top left and right panels, respectively), the bulge velocity dispersion and mass 
(central left and right panels, respectively), the cii'cular velocity of the galaxy (bottom left panel) and the virial mass of the DM halo (bottom right panel). Blue 
dots represent the outputs of the DeLucia2006a catalogue, grey and yellow shaded areas show the best fit to the model predictions and to th e observational 
datase ts, respectively. Starting from the u pper left panel down to th e lower right, the yellow shaded areas ref er to the best-fit relations obtained bv lMarconi et alJ 
i20Q4l) (the upper two panels of the plot). lFen-arese & For j i2003).lHaring & Rixl i2004l).lBaes et alji2003h and, in the lower-right panel, the four cui-ves show 
the equations 4 (cyan), 6 (green) and 7 f magenta) of lperraresd i2002l) and the results of Baesetal J i2003h (red). 



for the purpose of having /Edd(f) ^ decreasing function of time. 
Finally, from equations l ll4b . Ml\ and l |19| l, we have: 

Mbh(0 = MBH,peak + ^ [{l+Ctf - l] , (20) 

Where A = i^^, B = ^ + 1, C = (^)"" i. To derive 



equation ( I20t we set /Edd.O = 1 for continuity. We also impose 
/sdd = 10^"' as lower limit for the Eddington factor. 

Figure [T] shows the evolution of /Edd(f) (top panel), M^\i{t) 
(central panel) and Lbol(0 (bottom panel) for an illustrative case 
of a BH of Mbh = 10^ Mq accreting a mass Maccr = 5 x 10^ Mq, 
starting at z = 3, in the three prescriptions considered. The three 
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Relation Normalization (a) Slope (p) Scatter Scatterconcctcd 



log{MBH)-MK 


-4.37(0.24) 


-0.52(0.01) 


0.68 


0.53 


log(MBH) - Mb 


-0.61(0.17) 


-0.43(0.01) 


0.62 


0.53 


log(MBH)-log(Oc) 


-0.26(0.16) 


3.82(0.08) 


0.42 


0.28 


log(MBH)-log(Mb„lge) 


-2.39(0.19) 


0.96(0.02) 


0.58 


0.50 


log(MBH)-log(Vc) 


-1.61(0.18) 


4.05(0.09) 


0.45 




log(MBH)-log(MDM) 


-8.61(0.42) 


1.35(0.04) 


0.50 





Table 1. Parameters of the linear fits to the scaHng relations shown in Figure[2] A congelation of the form y = a + p • ,r has been assumed for all relations. The 
uncertainties in the normalizations and in the slopes are shown in parentheses. For details about the computation of the Scatter and the Scattercorrected see Sect. 

ED 



Relation 


a 


P 


Y 


Scatter 


l0g(MBH) -Mk 


17.29(0.10) 


1.25(0.01) 


0.04(0.01) 


0.51 


log(MBH) -Mk 


9.81(0.03) 


0.63(0.01) 


0.03(0.01) 


0.47 


log(MBH)-log(M|,„,„e) 


14.16(0.07) 


-2.21(0.01) 


0.15(0.01) 


0.44 



Table 2. Parameters of the fits to the scaling relations shown in Figure|3] A correlation of the form v = a + P •.v + y-.x" has been assumed for all three relations. 
The uncertainties in the parameters are shown in parentheses. For details about the computation of the Scatter see Sect.l3Jl 



green curves refer to lightcurve model ///, in which we set ^ = 0.5 
(short dashed), = 0.7 (dot-long dashed) and = 0.9 (short dashed- 
long dashed). 

Due to the present uncertainties concerning the origin of 
the BH seeds and their mass distribution, we assume MeH.seed = 
10"' M,;, for all seed BHs, irrespective of their halo host properties 
and their origin. Our results are robust with respect to this hypothe- 
sis since, as we have verified, they are basically unaffected by vary- 
ing AfBH.seed in the range [10^ — 10^]Mq at z < 3. More significant 
differences occur at higher redshifts, which we will investigate in 
detail in future work. 

The main parame t ers o f our model are listed in Ta- 
ble 1 of 'Croton et al.' | 2006|) , with the exception of, as in 
iDe LucisT ife Blaizot (2007), the values for the quiescent hot gas BH 
accretion rate, Kagn (defined in section I2.3.U . the star formation 
efficiency asF of equation and the instantaneous recycled frac- 
tion of star formati on to the cold disk, /?, which we set equal to 0.43 
(see Section 3.9 of lCroton et ai]|2006l) . 



3 MODELS VS. OBSERVATIONS 
3.1 The BH scaling relations 

Several observational evidences indicate that the masses of the BHs 
hosted at the centres of galaxies strongly correlate with different 
properties of their host bulges and DM haloes. In this section we 
compare the most recently observed BH scali ng relations at z = 
with t he predictions of the original model of IDe Lucia & BlaizotI 
( l2007l) . i.e. the predictions we obtain when assuming instantaneous 
mass accretion. We explore the effect of specifying the mass accre- 
tion rate at the end of this section. 



3.1.1 One-parameter relations 

In Figure |2] we show the correlation between the masses of the 
model BHs with six properties of their hosts, the K- and B-band 



bulge magnitude (Mb and M^), the bulge mass and velocity dis- 
persion (Mbuige and Oc), the circular velocity of the galaxy and the 
virial mass of the DM halo (V^ and Mdm). The blue dots represent 
the outputs of the model, while grey and yellow shaded areas show 
linear best fits to the model predictions and to the observational 
datasets, respectively. 

The dots in the plot refer to the population of BHs hosted in 
the central galaxies of FOF groups, or subhalos. We do not include 
those in satellite galaxies since in this case the host properties can- 
not be determined accurately. The data we have considered ar e: the 
Mbh - Mb and Mbh - Mk relations of 'Marconi et al.' ('20041) (top 
panels) the Mbh — f^c relation of Ferrarese & Ford (2005) (central 
left) the Mbh — Mbulge relation of Haring & Rix (2004) (central 
right) and the Mbh - Vc relation of ,Baes et al. (20_()3) (bottom left). 
No direct observational estimate is available for the Mbh — Mdm 
relation shown in the bottom right panel. The curves shown in 
this panel have been derived using different assumptions for the 
Mdm — Vc relation. In particular, the cyan, green and magenta lines 
correspond to equations (4 ), (6) and (7 ) of IFerraresd [2002) . while 
the red curve is taken from lBaes et al.l 1 20031) . 

Model predictions for Vc and Oc have been obtained by adopt- 
ing two different assumptions: i) Vc = V^ax, where Vmax is the 
maximum rotational velocity of the subhalo h osting the gal axy at 
its centre, and ii) Vc = l.SVvh- as derived bv ' SeliakI f2002h . The 
bu^ge velocity dispersion Oc is derived from the Vc — Oc relation of 
ISaes et al.l ( |2003|) . In the bottom panels, the grey areas correspond 
to a circular velocity obtained through hypothesis i) while the green 
ones, in better agreement with the data, assume hypothesis ii). 

The linear fit to the model data has been obtained using the 
bisector modification t o the ordinary least square s minimization 
approach, proposed by lAkritas & Bershadvl ( Il996h . for which the 
best-fit results correspond to the bisection of those obtained from 
minimizations in the vertical and horizontal directions. The estima- 
tor is robust and has the advantage of taking into account the pos- 
sible intrinsic scatter in the relation. The values of the best fit slope 
and the normalization are listed in Table [T] along with the scatter 
around the best fitting line. The uncertainties of the best fit param- 
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eters, also reported in the table, have been obtained by imposing 
Xd.o.f. = 1- 

As can be seen in Figure |2] the best fits to the model agree 
well with that to the data, within the scatter We note that, in all 
relations plotted, the scatter in the model is larger than that of the 
real data and also larger than the internal scatter observed in similar 
relations obtained from t he recent hydrodynam ical simulations of 
galaxy mergers (see e.g. Hopkins et al. 2007a). However, we no- 
tice that a large fraction of our model BHs are found in low-mass 
systems for which the scatter in the scaling relation is large. On 
the contrary, in the real datasets (and hydro-simulations) the major- 
ity of BHs belong to massive galaxies for which, according to our 
model, the scatter in the scaling relation is significantly smaller. To 
investigate whether the difference in the intrinsic scatter is real or 
is induced by a different sampling of the BH population, for each 
BH scaling relation we have discretized the range of the observed 
host galaxy properties in finite bins and generated 500 sub-samples 
by randomly extracting A'obs('^A:) model BHs from the parent cat- 
alogue, where Nobsi^x) is the number of BHs in the real dataset 
in each bin Ax . We have repeated the same fitting procedure in the 
500 sub-samples and found that the scatter is significantly reduced 
in this exercise, as indicated in the last column of Table[T] that lists 
the average scatter in the sub-catalogues. Therefore, the mismatch 
in the scatter results from sampling different BH populations: small 
objects in the model, massive objects in the observations. More- 
over, for the Mbh — c^c relatio n the scatter is very cl ose to 0.21, 
which is the value measured bv lHopkins eT^ j2007ah both in the 
observed and simulated data. 



3.1.2 Non-linear fits 

The agreement between model and data is satisfactory. However, 
we need to keep in mind that the model predictions for Vc and 
Oc and the observed relation between log(MBH) and log(MDM) 
have been obtained assuming further theoretical hypotheses. Con- 
sequently, the more constraining and reliable relations are the ones 
between the BH masses and the bulge magnitudes and masses. Fo- 
cusing on these relations and thanks to the huge number of model 
BHs, we have been able to investigate whether a non-linear fit 
provides a better match to the data. We find that the best fit is a 
quadratic function, y = a + ^ ■ x-\-y ■ . Figure [3] shows this fit 
(heavy green lines), together with the medians, the first and third 
quartiles (black points with error bars) of the model output, com- 
puted in a discrete number of bins. The internal scatter is signifi- 
cantly smaller than in the linear fit case. The values of the best fit 
parameters are reported in Table |2] While we predict, on average, 
too low BH masses for a fixed Mg with respect to the observa- 
tions (still consistent within the errors) the model predictions are 
in very good agreement with the data for the log (Mbh) and 
log(MBH) — log(Mbuige) relations. Interestingly, the 3-parameters 
fit of the la t ter rel ation is in excellent agreement with the one found 
bv lWvithel ( l200^ (magenta solid line in lower panel of Figure|3j. 



3.1.3 The fundamental plane relation 

In Figure|4]we compare the BH fundamental plane relation of our 
model at different redshifts with that obtained by iHopkins et al.l 
l l2007ah using both observational data and the outputs of hydro- 
dynamical simulations of galaxy mergers: 

log(MBH/M0) = 7.93 + 0.721og(Mji) + 1.41og(o2oo), 




-20 -22 -24 -26 
Mk (mag) 




108 lo'o lO'i 1012 



Figure 3. The tree model scaling relations best constrained by observa- 
tions. Here the black dots (with error bars) represent the medians (and the 
corresponding first and third quartiles) of the model outputs, computed in 
a discrete number of bins. The green lines show the best three-parameters 
fits to the model outputs (blue points). The magenta line in the lower panel 
refers to the best-fit relation obtained bv .Wvithe ( 2006.) . 



where Mjj is the galaxy stellar mass in units of 10 Mq, and 0200 
is the bulge velocity dispersion in units of 200 km s^^ . The red 
lines, bisecto rs of the plots , show the fundamental plane relation 
proposed by iHopkins et al.l ( l2007al) . Model prediction are repre- 
sented by blue dots, the black line is the best fit to the model and 
the shaded area its lo scatter. At low redshifts the agreement is very 
good. This is not surprising since at z ~ our model agrees with 
the Mbh ^ M^uige and Mbh ^ '^c scaling relations that represent 
fundamental plane projections. A discrepancy appears at high red- 
shifts. However, at z > 3 the fit involves only few objects and there- 
fore may not be very significant, especially when we account for 
the non-zero intrinsic scatter in the fundamental plane proposed by 
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Figure 4. The BH fundamental plane in the redshift range 0.1 ^ z ^ 5. The blue dot s are the model outputs, while the grey shaded areas show the best-fits to 
them. The red lines, con'esponding to the bisectors of the plots, are the predictions of lHopkins et alj J2007af) . The galaxy stellar mass, , is given in units of 
1O"M0, while the bulge velocity dispersion, 0200. is in units of 200 km 



iHopkins et al ] ( l2007ah . A remarkable success of our model is that 
it predicts very little evolution of the fundamental p l ane rel ation, 
at least out to z = 3, in agreement with lHopkins et all JlOOTah . The 
intrinsic scatt er, which does not evo lve with time either, is 3 times 
larger than in iHopkins et all ( l2007ah (we found a value around .6 
at all redshifts, while the one reported by IHopkins et al ] (l2007ah IS 
about 0.2). As discussed previously, the mismatch is reduced when 
using a number of model BHs consistent with the observed one. 



3.1.4 Dependence on the accretion history 



All scaling relations predicted by our model assume that BHs ac- 
crete mass instantaneously after merging events. What happens if 
we relax this assumption and specify the mass accretion rate in- 
stead? Figure |5] shows the impact of adopting different accretion 
recipes on the Mbh — Mbuice relation. As usual, filled dots repre- 
sent model predictions, grey shaded areas show the linear fit to the 
DeLucia2006a model scaling relation and the other hatched areas 
indicate the linear fit to the model predictions obtained with our 
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Figure 5. The log(MBH) — log(Mbuige) scaling relation for our different pre- 
scriptions for the BH mass accretion. The filled dots represent model predic- 
tions, the grey shaded areas show the linear fit to the DeLucia2006a model 
scaling relation and the other hatched areas indicate the linear fit to the /, 
// and /// lightcurve models, as indicated by the labels. The black dots and 
grey shaded areas, in the lower right panel, show the prediction obtained 
with the parameterization given by the equations 12 It . as explained in Sec- 
tion[331 

different recipes, as indicated by the label Clearly, these predic- 
tions depend little on the assumed mass accretion histories for each 
individual quasar event (the fit parameters have fluctuations of no 
more than about 1%). This is a consequence of the fact that the 
BH scaling relations depend mainly on the total mass accreted, and 
very little on the time spent in the accretion process. We have veri- 
fied that all other scaling relations, including also the fundamental 
plane relation, does not change significantly by adopting any of the 
mass accretion prescriptions described in Section [2.3.3l 

3.2 The BH mass function 

The BH mass function (MF) is defined as the differential co- 
moving number density of BHs as a function of their mass. In 
Figure [6l we compare the BH MF predicted by our model for the 
prescripti ons / (blue line), // ( red) and /// (green) with those ob- 
served bv lShankar et alj |2004|) (grey area) and by Shankar (2007, 
in preparation) (yellow ar ea) at z ~ 0. In neither case the BH masses 
were determined directlvi Fshankar et al. |(|2004^ derive the BH mass 
from the observed Mbh — ^hiil ge relation v vhile S hankar (2007) 
use the Mbh — f^c relation of iTundo et alj ( l2007h . We note that 
the model BH MF is in good agreement with the observed ones, 
within the mass range accessible to observations exept in the inter- 
val ~ 10^ — lO^Mf), in which the number density of model BHs is 
smaller than the observed one. 

The reason of the small mismatch between the observed and 

' The meaning of the black dots and shaded areas in the bottom left panel 
of Fig. 5 is discussed in Section [33l 
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Figure 6. Comparison of the BH mass function predicted by lightcurve 
model s /, // and /// with the one observationally derived by IShankar et alJ 
(20041), and with the new one obtained by Shank ar (2007, in preparation) 
using the Mbh — o relation bv lTundo et alj j2007h . The grey areas show the 
prediction obtained with the parameterization given by the equations i2l\ . 
as explained in Section [331 



the model BH MFs will be investigated in a forthcoming paper in 
which we study the redshift evolution of the BH MF and its depen- 
dency on the properties of the host galaxy. Finally, we note that, as 
shown in Figure [S] the model predictions for the BH MF are ro- 
bust with respect to the prescription adopted for the mass accretion 
history of the individual quasar episodes. 



3.3 The AGN bolometric luminosity function 

The luminosity function (LF) of AGN, namely the derivative of 
their co-moving number density with respect to luminosity, rep- 
resents a unique tool to understand their cosmological evolution. 
Semi-analytic models predict the total (bolometric) luminosity of 
a statistically complete AGN catalogue, and to compare model 
LFs with observations we need to specify a bolometric correc- 
tion, i.e. how to convert th e luminosities ob s erved in a particular 
band into bolometric ones telvis et alJ[T994l : iMarconi et id]|2004 
'Hop kins etai] |2006^. Another correction is required to acc ount for 
possi ble incompleteness effects (see e.g. Comastri 2004; Gil li et all 
l2007i) , which includes the possible existence of a population of ob- 
scured AGN whose fraction may depend on the wavelength band 
and redshift (lElvis et al.ll994lMarconi et al.l2004lLa Franca et al.l 
l2005l;lLamastra et al.ll2006h . 

Here w e compare our predic tions with the bolometric LF 
obtained by iHopkins et alj ( l2007h from the LFs observed in 
diffe rent bands: radio (see e.g. iNagar et al |2005|), optical (see 
e^ 



1997 



Kennefick et al 



iGrazian et al 



Schmidt et alT I995I : iKoehler et al.l 



'Fanetar'~'"2001 



Hunt et alj 120Q4 
Richards et al. 2005, 2006; Siana 
iShankar & Mathun i2007i; jBong iorno 



'Wolf et al' '2003|; 



Cristiani et al. 2004; Groom et al. 2005 



et a kbooej; iFontanot et al.H2007 
srno et al.' '2007*), infra-red (sei 



e.g. iBrown et alj|2006l: iMatute et al. 2006; Babbed ge et al.ll2006h. 
soft X-ray (see e.g. iMivaii et al.l I200QL 120011 ; [Silverman et al.l 
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Figure 7. The bolometric LFs predicted by our lightcurve models / (blue bands), // (red bands) and /// (green bands), in the redshift range 0.1 5C ^ 5. are 
here compared with the best-fits to observational data obtained bv lHopkins eTaU feo07) (yellow bands). The grey areas show the predictions obtained with 
the parameterization given by the equations j2U . as explained in Section [33] Uncertainties in the model LFs are computed by assuming Poisson statistics. 
The dashed vertical gre en lines mark the ran ge of the bolometric luminosities accessible to observations. The dotted red vertical lines show the luminosities 
beyond which the LP of lHopkins et al predicts a number of AGN in the whole volume of our simulation smaller than 1 0. The vertical grey dotted lines 

around the red ones have been calculated considering the en'or in the best-fit of lHopkins et al.l i2007l) . 



2005 




Hasineer et alj 2005b, hard X-ray (see e.g. Bareer et al. 


2003 


UedaetalJ 20031: iBarser et al.l l2003l: iNandraetal. 


2005 


Sazonov & RevnivtsevI |2004 Isilverman et al.l l2005d; 


La Franca et al. 20051: Shinozaki et al.l |2006|; Beckmann et al. 


20061). and from emission lines (see e.g.lHao et alj|2005'). Uncer- 



tainties in these corrections contribute to the scatter in the observed 
LF, i.e. to the width of the yellow areas in F igure [7] that show the 
AGN bolometric LF of lHopkins et al.l (120071) at different redshifts. 



The model predictions are also represented by areas with different 
colours, with a width corresponding to lo Poisson error bars. 
The vertical, green dashed lines bracket the bolometric luminosity 
range accessible to observations. The vertical, re d dotted lines 
show the luminosities beyond which the LF of Hopkins et al.l 
(^007) predicts less than 10 AGN in the volume of our simulation, 
i.e. the maximum luminosities at which our model BH sample 
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is statistically meaningful; lo uncertainties on this maximum 
luminosity are represented by the two grey dotted lines. 

From Figure [7] we see that, on average, type-/ lightcurve un- 
derestimates the AGN number density at all epochs. However, 
while at high redshifts the model matches the faint-end of the LF 
and underpredicts the number density of the bright objects, the sit- 
uation is completely reversed at z ~ 0, where the model correctly 
reproduces the number density of bright AGN but underestimates 
the faint ones. At low redshifts the problem can be alleviated by re- 
ducing the Eddington factor, as in our type-// lightcurve. However, 
in this case the discrepancy between model and data at high red- 
shifts increases. Adopting the type-/// lightcurve allows to match 
observations in the whole range of luminosities in the redshift range 
0-5 ^ z 1> but overestimates the number of luminous AGN at 
z < 0.5 and underestimates them at z > 1. 

Therefore, we conclude that in our present semi-analytical 
framework we can reproduce the observed AGN LF at low and in- 
termediate redshifts. However, at z > 1, we under-predict the num- 
ber density of bright AGN, regardless of the BH mass accretion 
rate and light curve model assumed for each quasar episode. To in- 
vestigate if it is possible to modify our prescription for the mass 
accretion to fit the AGN LF at all redshifts, we tried different val- 
ues of /Edd and e as a function of / and Mbh, within physically 
motivated ranges. Despite of the considerable freedom in choosing 
fEAiiU^^u) we failed to find a model able to match simultane- 
ously the observed BH scaling relations, the BH MF, and the AGN 
LF, especially at high redshifts. We also used different plausible 
values for the BH seed mass, and we still were not able to fit the 
high-z LF. We interpret this failure as an indication that our theo- 
retical framework itself is inadequate to account fully successfully 
for the AGN phenomenon. 

One possible way out is to modify the model assumptions for 
the efficiency of BH growth in the quasar mode following mergers 
at high z. A significant improvement of our results at high redshifts 
can for example be obtained by substituting equation M2\ and ( I13t 
with 



f /bh = 0.01 ■ log (jf^ + l).z z> 1.5andMBH>lO''M0 
\ AMbh.q = 0.01 mcoid z>6 

(21) 

while keeping prescription /// for the quasar lightcurves. The pre- 
dictions of this new model for the log(MBH) — ^og(M\,^\„^) scaling 
relation is shown as black dots in the bottom-right plot of Fig. [5] 
Model predictions for the BH MF and AGN LF are shown in Fig- 
ures [6] and [T] respectively. An accretion efficiency that increases 
with the redshif t has been already advocated in the dynamical 
model of Croton (2006). A physical justification to this assumption 
is provided bv iMo et ah ( il998i) . Indeed, their model predicts that 
galactic disks were more centrally concentrated in the past, making 
it more efficient the BH feeding at high redshift. It is worth stress- 
ing that equation J21l ( might not provide the best fit to the data as 
we did not explore the parameter space systematically. However, it 
suggests that a good match to the observed scaling relations, BH 
MF and AGN LFs can be obtained within our semi-analytic frame- 
work by modest changes of the BH growth at high redshifts. The 
solution provided by equation J2U is not unique either, since larger 
amounts of mass can be accreted also by invoking alternative mech- 
anisms that trigger gas accretion episodes, for example by secular 
evolution through disk ins tabilities, or by al luding to a higher gas 
cooling efficiency (see e.g. lViola et alj|2007l) . 



4 CONCLUSIONS 

In this paper we have used and extended a semi-analytic model for 
the co-evolution of galaxies and their ce ntral BHs, developed on the 
outputs of the Millennium Simulation (Springel et al. 2005), and 
descr ibed in detail in Croton et al. (2006) and De Lucia & BlaizoJ 
( |2007|) . The aim of the model is to reproduce the observed proper- 
ties of BHs, AGN and their galaxy hosts. The physical assumptions 
in the model with respect to BH growth can be divided into two 
sets. The first one concerns the mass accretion history of the cen- 
tral BHs in h alos, where we disti nguish between radio mode and 
quasar mode jCroton et alj|200^ . This set makes predictions for 
the relation between BH and galaxy host properties, which can be 
compared to the observed scaling relations between BH mass and 
different properties of their host galaxies. The second set of pre- 
scriptions specifies the detailed AGN activity and lightcurve of in- 
dividual quasar episodes, and leads to predictions for the AGN LF 
as a function of redshift. We considered three different models for 
this detailed AGN activity, one of them motivated by the results of 
recent hydrodynamical sim ulations of galaxy me r gers that include 
BH growth and feedback jHopkins et^ bOOSlJpi Matteo et al.1 
120051 : ISpringel et alj2005l) . 
The main results of our analysis are as follows: 

(;) The semi-analytic model is approximately able to repro- 
duce the observed BH scaling relations over the whole range of 
BH masses and galaxy properties probed by observations. The in- 
trinsic scatter in the model is significantly larger than in the data, 
a mismatch that can in part be accounted for by adopting the ob- 
servational selection criteria to obtain a mock BH catalogue with 
similar characteristics as the observed one. 

{ ii) We find evidence that a quadratic relationship provides a 
significantly better fit to some of the model scaling relationships 
than a linear one, as already noticed bv lWvith^ ( l2006i) . 

(Hi) Our m odel also matches the BH fundamental plane re- 
lation derived bv lHopkins et al.l ( 12007 al) . and successfully predicts 
very little evolution of this plane, at least out to z ~ 3. 

(iv) The model BH mass function is in good agreement with 
the observed one within the mass range accessible by observations, 
except on the range ~ 10^ — 10^ Mp,, in which the number density 
predicted by the model is smaller than the observed one. 

(v) Model predictions for the BH mass function, scaling rela- 
tions and fundamental plane relation are basically unaffected when 
using different prescriptions for the AGN lightcurves of individual 
quasar events. This is because these predictions are only sensitive 
to the model assumptions for the absolute growth of the BHs in 
each merger event. 

(vi) The AGN LF is systematically underestimated by as- 
suming that BHs accrete mass with a constant Eddington factor 
/Edd = 1 • The detail of the discrepancy, however, change with red- 
shift since at high z the model matches the faint-end of the LF but 
underpredicts the number density of the brightest objects, while 
the situation is reversed at z ~ 0, in agreement with the results of 
several semi-analytic models (see, e.g. Marulli et al. 2007, and ref- 
erences therein). Reducing the Eddington ratio, as in our lightcurve 
model //, alleviates the faint-end mismatch but amplifies the bright- 
end discrepancy at high redshifts. A significant improvement at low 
redshifts is obtained when the Eddington-limited growth of the BH 
is followed by a lo ng quiescent phase vy ith lower Eddington ra- 
tios, as suggested bv lHopkins etalj|2005l) and implemented in our 
lightcurve model ///. In this case our model is able to match the ob- 
served AGN LF in the interval 0. 1 < z < 1, over the whole range of 
luminosities that are accessible to observations and where our pre- 
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dictions are statistically significant. However, our predicted number 
density of bright AGN is still biased low at z > 1 . 

(vii) Our model is able to account for all observations consid- 
ered in this work apart for the AGN LF at high redshifts. We were 
not able to eliminate this mismatch by simply modifying the accre- 
tion efficiency, e, the Eddington factor, /Edd^ or the BH seed mass 
(when considered in physically plausible ranges). Clearly, we need 
to modify assumptions in the underlying semi-analytic framework 
for BH growth. A simple, ad hoc increase of the mass fraction ac- 
creted during the quasar mode at high redshift can indeed remedy 
the problem. However, this solution is not unique as several high- 
redshift modifications to the original model, like new mechanisms 
that trigger BH activity in addition to galaxy merging or more effi- 
cient gas cooling resulting in a larger reservoir of cold gas, can be 
advocated to bring the predictions in line with observations. How- 
ever, it remains to be seen whether any of these alternatives is phys- 
ically plausible. 

f viii) Our model predictions at z < 3 are robust to changes in 
the assumed BH seed mass, but are sensitive to it at larger redshift. 
We will further explore this issue in a subsequent paper where we 
plan to study to what extent current observations can constrain the 
seed BH mass function. 

From our analysis we conclude that the AGN LF at high red- 
shifts constitutes a strong constraint for semi-analytic models that 
describe the co-evolution of galaxies, BHs and AGN. This sug- 
gests that significant improvements can be obtained in two ways. 
From the theoretical side, we need to develop a physically mo- 
tivated mechanism capable of increasing the number density of 
bright AGN at z > 1 without modifying the model predictions at 
low redshifts. From the observational point of view, we need to 
improve the AGN LF estimates at high redshift, both by enlarging 
current high-z AGN samples and by reducing the current uncertain- 
ties originating from bolometric and incompleteness corrections, in 
particular for the population of Compton Thick AGN. In addition, 
other observational tests should be performed, like the ability of 
our model to match the observed AGN clustering, as quantified by 
the angular and spati al two-point correlations function. In particu- 
lar, iLidz et al] ( l2006h pointed out that the luminosity dependence 
of quasar clustering can discriminate between different lightcurve 
models, a question we will address in a forthcoming work. 
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